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Since its introduction in 1976, 1 the Rys Quadrature method has proven a very 
attractive method for evaluating electron repulsion integrals for calculations using 
Gaussian type orbitals. Since then, there have been considerable refinements of the 
method, 2 ’ 3 but at it’s core, Gaussian weights and nodes are used to exactly evaluate 
using a numerical approach the transform integral. One of the powers of the Rys 
Quadrature method is the relative ease in evaluating integrals involving functions of 
high angular momentum, and Seward 2 can handle up to i functions (/ = 6). This 
requires Rys quadrature of order 2-H + 1 = 13, yet when we tried to compute 
weights and nodes for higher order quadratures, we ran into numerical difficulties. 
In this work we report on the complete resolution of these numerical difficulties, 
and we have easily computed accurate quadrature weights and nodes up to order 
101, thus now functions up to F = 50 can be handled, if desired. All calculations 
were carried out using 128-bit precision. 

King and Dupuis very nicely described the calculation of Rys quadrature weights 
and nodes, 4 and we initially followed many of their procedures. However we used 
more sophisticated algorithms for computing the weights and nodes. In our first 
attempt, we computed the moments 


for the fixed parameter X in terms of the incomplete gamma function via the 
algorithm given by Press et al. s and then used them with the algorithm of Wheeler, 6 
an efficient reformulation of the Sack and Donovan 7 method, to compute the weights 
and nodes. In order to compute a jY point quadrature, the moments for 

a = 0,1 2jY — 1 are required. However, even though we used 128-bit precision, 

we had difficulties in obtaining all positive nodes for the higher order quadratures. 

In our second attempt, we computed the modified moments 
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where ;t£ is an orthogonal polynomial. Note that now the integration interval has 
changed to be symmetric about the origin - this means that all odd ii lead to zero 
modified moments, so for a jY point quadrature, the moments for a = 0,1, ... ,4jY are 
required. In the limit X 0, the Rys polynomials become Legendre polynomials, and 
so for "small” X we used a=L, i.e., the jt£ are Legendre polynomials. In the limit 
X oo, the Rys polynomials are proportional to Hermite polynomials of argument 
t/\'X, and so for "large” X we used a=H, i.e., the jt£ are Hermite polynomials of 
argument t/\X. By means of these choices, rather than increasing in magnitude as n 
increases, the modified moments rapidly decrease as n increases. But the trade off is 
the modified moments must be evaluated numerically. Initially we tried Romberg 
integration, 5 whereby the integral was evaluated by the trapezoidal rule for step 

. n. n. h. 

sizes n h - 1 - , - , ... and then extrapolated to zero step size. This eliminated the non- 
positive node problem seen in the first attempt, and we moved onto the next stage 
of the development, namely the interpolation of the weights and nodes with respect 
to*. 

We follow King and Dupuis 4 and interpolate using a Chebyshev expansion of the 
dependence of the weights and nodes on X. This turned out to be a critical factor in 
obtaining accurate high order weights and nodes. This is because the Chebyshev 
expansion coefficients are a very good indicator of numerical noise in the 
calculations, and we found that although our weights and nodes looked reasonable 
for individual X, they clearly showed unacceptable noise levels when looked at as a 
function of X. This lead to our next level of refinement. We experimented with 
several ways of carrying out the Chebyshev expansion and ended up dividing X into 
evenly space intervals of width unity. 

In order to reduce the numerical noise, we tried increasing the accuracy of the 
computation of the modified moments, but found that the Romberg integration 
simply could not provide the required accuracy. The first step in improving the 
process was to switch from the trapezoidal rule to repeated 10 point Gauss- 
Legendre quadrature. We evaluated the integrals by dividing the range into 
l,Z r i r S r equal segments, integrating each segment by the 10 point quadrature rule, 
then stopped when two successive number of segments gave the same results to 
within an input error tolerance, typically 10 This turned out to be much more 
efficient than the Romberg integration, but still was not able to deliver the required 
accuracy. The solution to this problem was to look at the integrand: the reason the 
desired accuracy was hard to obtain was that the integrand was highly oscillatory 
and thus there was significant cancellation between various regions. Thus, we then 
determined the nodes in the integrand and split up the integration range so that 
each range started and ended on a node, and included a node between the starting 
and ending range, if possible. This procedure worked very well. We determined the 


nodes via bisection, 5 with the bisection search stopping when the node was 
determined to nearly machine precision (10 -:u ). 

By means of this procedure, we could generate accurate weights and nodes for 
quadrature’s up to jV = '18 before unacceptable numerical noise again crept into the 
calculations. This time we attribute the problem to the fact that the modified 
moments for high order become so small that numerically they were similar to the 
odd modified moments that are all identically zero. 

This then led to the final procedure. If one reads the discussion of algorithms in Ref. 
5, a primary driver is numerical efficiency, with numerical stability being an 
important, but secondary concern. In the present application, numerical efficiency is 
not an important concern, since the quadrature weights will be computed only once. 
Thus, we will proceed as follows: we will directly determine the coefficients and 
i? n in the recursion relation 

ifr*) = (t- - t*) - Mu-iCrO 
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Eq 4 

for ii -2 = n and n — 1 subject to the boundary conditions ■_ =0 and ii u = 1, where 
we imply the dependence of a n and i? n on X. This is the procedure of Stieltjes, 5 and 
requires the evaluation of the integrals and 

(Rj,-! J. These are evaluated using the repeated 10 point quadrature described 
above. Once the recursion coefficients are determined, we determine the nodes by 
root finding, using the intervals from the roots of the previous order polynomials as 
initial guesses to a several step bisection search followed by Newton-Raphson 
refinement. Since the functions are evaluated by recursion, the derivative required 
for the Newton-Raphson procedure is easily computed and thus we have all the 
required information for evaluating the quadrature weights via 5 

tfth-iilft-i), 
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Eq 5 


where ft is one of the roots of JJ lV . 

This method requires about 'ii V integrals for a jY point quadrature, whereas the 
previous algorithm requires about ZjV integrals, so is clearly less efficient, but it 
gives very good results, for the recursion coefficients are given as a simple ratio of 
integrals: a„ = andi^ = (Rn-i^n-i)- We have used 


this method up to jV = 101 with absolutely no sign of numerical noise at the 10 -it 
level. 

It is instructive to compare the present method with that used by King and Dupuis. 1 2 3 4 5 6 7 
The principle difference is they based their calculations on the moments of Eq 1, and 
used Schmidt orthogonalization to determine the Rys polynomials. Our procedure, 
in contrast, explicitly evaluated the integrals and (R,,^ J at 

each step once R 1; had been determined. This is exactly the difference between the 
Gram-Schmidt orthogonalization and Modified Gram-Schmidt orthogonalization, 8 
and is the reason Modified Graham-Schmidt orthogonalization has such superior 
numerical stability compared to the original version. 

We now address the issue of evaluating the quadrature weights and nodes in the 
large X limit. As mentioned previously, as X oo, the Rys polynomials become 
scaled Hermite polynomials, thus the Rys quadrature weights and nodes can be 
computed from Gauss-Hermite quadrature weights and nodes. King and Dupuis 4 
examined this transition carefully, and came up with a fitting procedure which 
enabled them to interpolate between the asymptotic limit and the X values used in 
the Chebyshev expansion. We have found that this smooth transition was present 
for low weights and nodes, v.\ and ft for i « jV, but not for high weights and nodes, 
Wi and ft for i ^ jV. Thus for simplicity, we carried out the Chebyshev expansion 
until the asymptotic region where the Rys quadrature weights and nodes can be 
computed accurately from Gauss-Hermite quadrature weights and nodes. Compared 
to the algorithm of King and Dupuis, 4 our method requires more memory, but this is 
not an issue for modern computers. 

We have every expectation that the developments reported herein will enable some 
very exciting new calculations in molecular electronic structure theory. 
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